Data Replication Analysis
Victoria Zdanowicz
11/24/2021
Assignment Overview
You do not need to replicate ALL of the analyses presented in the paper, but at minimum you must replicate at least 3 analyses, including at least one descriptive statistical analysis and one inferential statistical analysis. As part of this assignment, you must also replicate to the best of your abilities at least one figure.
Stargazer package -> good tables, publication ready kable -> knit into HTML DT -> allows you to navigate table in RMarkdown
library(curl)## Using libcurl 7.64.1 with LibreSSL/2.8.3library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(broom)
library(ggplot2)
library(gridExtra)##
## Attaching package: 'gridExtra'## The following object is masked from 'package:dplyr':
##
## combinelibrary(sciplot)
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodelibrary(glmmTMB)
library(DHARMa)## This is DHARMa 0.4.4. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')library(emmeans)
library(sjPlot)
library(huxtable)##
## Attaching package: 'huxtable'## The following object is masked from 'package:sjPlot':
##
## font_size## The following object is masked from 'package:ggplot2':
##
## theme_grey## The following object is masked from 'package:dplyr':
##
## add_rownamesI first selected an article titled, “Attractiveness of female sexual signaling predicts differences in female grouping patterns between bonobos and chimpanzees” by Surbeck et al. (2021) to replicate for this assignment. I found the study and their results extremely interesting, especially since I recently heard Surbeck present about his research at NEEP and have long heard the idea that bonobos are more gregarious than chimps because of ecological differences. However, as I read through the paper and began to try and manipulate the dataset they shared on Dryad, I realized I was missing a lot of the details necessary to run the models they did in their analyses. The results and methodology sections at the end of the paper fail to include the parameters of their GLMMs, making it nearly impossible for me to try and understand and replicate - especially because this kind of modeling is brand new to me!! I then decided to look back over the other articles I had been considering for the assignment and realized the Lee et al. (2021) paper titled, “Gregariousness, foraging effort, and affiliative interactions in lactating bonobos and chimpanzees” covered a somewhat similar topic as the Surbeck paper while including much more detailed descriptions of their analyses!
Getting Started
First loaded in the data. The dataset on Dryad was saved as an Excel file with 2 sheets: one with individual chimp/bonobo data on time spent ‘feeding, traveling, social’ and the other with individual time spent ‘alone.’ I exported each individual sheet into separate .csv files and uploaded them to GitHub so that I could curl the data into R.
#loading in 'alone' data
a <- curl("https://raw.githubusercontent.com/vzdanowicz/AN588_DataReplication_vrz/main/lee_alone.csv")
time_alone <- read.csv(a, header = TRUE, sep = ",", stringsAsFactors = FALSE)
names(time_alone)<-c('species','ID','age','total','alone')
#loading in 'feed/travel/social' data
b <- curl("https://raw.githubusercontent.com/vzdanowicz/AN588_DataReplication_vrz/main/lee_behav.csv")
feed_travel_social <- read.csv(b, header = TRUE, sep = ",", stringsAsFactors = FALSE)
names(feed_travel_social)<-c('species','ID','age','total','feed','travel','social','social_adj')
#view
head(time_alone)| species | ID | age | total | alone |
|---|---|---|---|---|
| bonobo | Gwe | 0--18 | 68 | 0 |
| bonobo | Iri | 18--36 | 88 | 0 |
| bonobo | Iri | 36--54 | 40 | 1 |
| bonobo | Nin | 18--36 | 110 | 0 |
| bonobo | Olg | 0--18 | 70 | 0 |
| bonobo | Olg | 36--54 | 84 | 1 |
head(feed_travel_social)| species | ID | age | total | feed | travel | social | social_adj |
|---|---|---|---|---|---|---|---|
| bonobo | Dju | 0--18 | 1618 | 662 | 316 | 272 | 137 |
| bonobo | Gwe | 0--18 | 3493 | 1415 | 386 | 720 | 334 |
| bonobo | Iri | 0--18 | 2794 | 993 | 471 | 579 | 307 |
| bonobo | Nin | 0--18 | 2621 | 1089 | 370 | 367 | 52 |
| bonobo | Olg | 0--18 | 1425 | 383 | 202 | 95 | 5 |
| bonobo | Uma | 0--18 | 2122 | 707 | 400 | 367 | 155 |
The article contains a fairly detailed description of their statistical analyses in R. I used the help() function to learn more about the glmmTMB package and function used by these authors. glmmTMB = fit generalized linear mixed models using Template Model Builder (TMB), formula follows lme4 syntax
formula,data = NULL,family = gaussian(),ziformula = ~0,dispformula = ~1,weights = NULL,offset = NULL,contrasts = NULL,na.action, se = TRUE,verbose = FALSE,doFit = TRUE,control = glmmTMBControl(),REML = FALSE,start = NULL,map = NULL,sparseX = NULL
Exploring & visualizing the datasets..
par(mfrow = c(1,2))
boxplot(data = time_alone, alone ~ age * species, group = time_alone$age)
boxplot(data = feed_travel_social, feed ~ age * species, col = c('cadetblue1','darkseagreen1'))
**do I need to convert species into factor?? in looking at data from m
#convert column 'id' from character to factor
time_alone$ID <- as.factor(time_alone$ID)
time_alone$species <- as.factor(time_alone$species)
time_alone$age <- as.factor(time_alone$age)
str(time_alone)## 'data.frame': 30 obs. of 5 variables:
## $ species: Factor w/ 2 levels "bonobo","chimpanzee": 1 1 1 1 1 1 1 1 1 1 ...
## $ ID : Factor w/ 19 levels "BAH","DL","EZA",..: 9 10 10 11 12 12 13 13 14 15 ...
## $ age : Factor w/ 3 levels "0--18","18--36",..: 1 2 3 2 1 3 1 2 1 3 ...
## $ total : int 68 88 40 110 70 84 48 27 30 119 ...
## $ alone : int 0 0 1 0 0 1 2 2 0 0 ...To test our three predictions (described above), we fitted generalized linear mixed models (GLMMs) to each response variable (response variables for each prediction described below) using the glmmTMB function in the glmmTMB package with a beta-binomial error structure. We initially fitted GLMMs using binomial error structures but found that all models were overdispersed. Overdispersion occurs when variance is higher than predicted by the model because the model lacks an adjustable dispersion parameter (e.g., as in binomial and Poisson models; Bolker et al. 2009; Zuur et al. 2009). Beta-binomial models include an adjustable dispersion parameter that allows the model to predict variance appropriately for binomial proportion data (Harrison 2015).
They reported results of nonparametric dispersion tests for all models using the testDispersion function (case sensitive) in DHARMa. None of the beta-binomial models exhibited overdispersion.
Model assumptions were evaluated by ’visually assessing QQ plots and the distribution of residuals plotted against fitted values using the simulateResiduals (case sensitive) function in DHARMa
Prediction 1: Time Alone
P1: Lactating female chimpanzees spend more time alone than lactating bonobos
Lee et al. (2021) shares the response variable (prop_alone) was calculated by dividing # alone_scans of a given female by the total # of party scans collected on that female during that age class. Thankfully, the data is already organized by each species, female ID, and infant age class, so I only had to manipulate my dataframe time_alone to calculate the necessary response variable.
To do this - I created a new column in time_alone and then wrote a for loop that filled the appropriate proportion values into this new column
`time_alone$prop_alone <- NULL
for (i in 1:nrow(time_alone)) { time_alone\(prop_alone[i] <- time_alone\)alone_scans[i]/time_alone$total_scans[i] }`
Now that I have my response variable, I can begin trying to work through the GLMM. I started this process by carefully reading through the article to understand the parameters used, as well as read the {glmmTMB} package information and ecological GLMM examples from our course supplementary readings. I also had to Google what exactly a beta-binomial model meant (as described in the source article).
I first tested the interaction between species and infant age class (as outlined in the article). I first triede to create a full model that included the interaction and then a reduced model with both as independent fixed-effects and test the relationship with Anova.. like we had done in Mixed Effects module. But I kept getting the error message Error in fixef(mod)[[component]] : invalid subscript type ‘list’
I then realized Anova() function from {car} can actually directly test interactions within a single model. So I paired down to just have P1_interaction model and then tested the signifiance of the interaction using Anova() Wald “Chisq” with type “III”.
P1 - Interaction (M1)
#testing species/age interaction
P1_interaction <- glmmTMB(prop_alone ~ as.factor(species) : as.factor(age_bin) + (1 | id), data = time_alone, family = betabinomial) Anova(P1_interaction, type = c("III"), test.statistic = "Chisq")
this includes individual effects ^^ but they only report interaction effects?
model baseline bonobo 0-18. can i run the model without creating new column?? this was much more successful… ALSO needed to add weight = total_scans… so weight with denominator of the proportion of the response variable
I could not figure out why my results did not look like those in the article… the values of the model above were close but didnt seem to be giving thee right output. I continued to read about the glmmTMB function and package to try and understand what piece I was missing. I eventually figured out that in order to have a proportion (non binary/Bernoulli) response variable in the model, it needs to be weighted. The {glmmTMB} package says: “Binomial models with more than one trial (i.e., not binary/Bernoulli) can either be specified in the form prob ~…,weights = N, or in the more typical two-column matrix cbind(successes,failures)~… form” (Magnusson et al. 2021).
I found an example of someone doing this in a question thread online. They kept the response variable as the proportion without creating an additional column in their data frame and used the denominator as the weight. After attempting that version of the model, my output was finally identical to Lee et al. (2021)!!
P1M1 <- glmmTMB(alone/total ~ species * age + (1 | ID), data = time_alone, family = betabinomial(link = "logit"), weight = total)
P1M1_anova <- Anova(P1M1, type = c("III"), test.statistic = "Chisq")look at model/anova results and compare with paper
summary(P1M1)## Family: betabinomial ( logit )
## Formula: alone/total ~ species * age + (1 | ID)
## Data: time_alone
## Weights: total
##
## AIC BIC logLik deviance df.resid
## 151.8 163.0 -67.9 135.8 22
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## ID (Intercept) 1.009e-08 0.0001005
## Number of obs: 30, groups: ID, 19
##
## Dispersion parameter for betabinomial family (): 13.6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.9000 0.7189 -5.425 5.8e-08 ***
## specieschimpanzee 2.5870 0.7646 3.383 0.000716 ***
## age18--36 -0.8451 1.2221 -0.692 0.489236
## age36--54 0.5001 0.9182 0.545 0.585949
## specieschimpanzee:age18--36 0.6804 1.3018 0.523 0.601210
## specieschimpanzee:age36--54 -0.8054 1.0334 -0.779 0.435795
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print(P1M1_anova)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: alone/total
## Chisq Df Pr(>Chisq)
## (Intercept) 29.4281 1 5.803e-08 ***
## species 11.4475 1 0.0007159 ***
## age 1.3684 2 0.5044876
## species:age 1.5098 2 0.4700581
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1need to pull out the chisq x2, df, and p val for species:age_bin … using {broom} (from glm module)
P1M1_results <- tidy(P1M1_anova)
P1M1_x2 <- P1M1_results$statistic[4]
P1M1_df <- P1M1_results$df[4]
P1M1_p <- P1M1_results$p.value[4]P1 - Independent (M2)
Interaction was not significant so we need to model them as independent fixed effects.. now testing independent effects so use type II anova instead of 3
P1M2 <- glmmTMB(alone/total ~ species + age + (1 | ID), data = time_alone, family = betabinomial(link = "logit"), weight = total)
P1M2_anova <- Anova(P1M2, type = c("II"), test.statistic = "Chisq")looking at model results to compare to article
summary(P1M2)## Family: betabinomial ( logit )
## Formula: alone/total ~ species + age + (1 | ID)
## Data: time_alone
## Weights: total
##
## AIC BIC logLik deviance df.resid
## 149.4 157.8 -68.7 137.4 24
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## ID (Intercept) 4.269e-08 0.0002066
## Number of obs: 30, groups: ID, 19
##
## Dispersion parameter for betabinomial family (): 13.9
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.8044 0.5085 -7.481 7.37e-14 ***
## specieschimpanzee 2.4699 0.4814 5.130 2.89e-07 ***
## age18--36 -0.2671 0.4194 -0.637 0.524
## age36--54 -0.1396 0.4048 -0.345 0.730
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print(P1M2_anova)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: alone/total
## Chisq Df Pr(>Chisq)
## species 26.3209 1 2.891e-07 ***
## age 0.4135 2 0.8132
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1need to pull out the chisq x2, df, and p val for species:age_bin … using {broom}
P1M2_results <- tidy(P1M2_anova)
#pulling out values for species
P1M2_SPECIES_x2 <- P1M2_results$statistic[1]
P1M2_SPECIES_df <- P1M2_results$df[1]
P1M2_SPECIES_p <- P1M2_results$p.value[1]
#pulling out values for age
P1M2_AGE_x2 <- P1M2_results$statistic[2]
P1M2_AGE_df <- P1M2_results$df[2]
P1M2_AGE_p <- P1M2_results$p.value[2]not 100% sure how yet to pull out variables from the model or anova.. need to assess what I need to create the table (tables 3/4 in the paper)
P1 - Dispersion Tests
‘We reported results of nonparametric dispersion tests for all models using the testDispersion function (case sensitive) in the DHARMa package. None of our beta-binomial models exhibited overdispersion. We evaluated model assumptions by visually assessing quantile–quantile plots and the distribution of residuals plotted against fitted values using the simulateResiduals (case sensitive) function in the DHARMa package.’
‘The nonparametric dispersion tests were not significant for either Time Alone model (interaction effect model: deviance ratio = 0.957, P = 0.960; independent effects model: deviance ratio = 1.002, P = 0.928).’
Next need to look at nonparametric dispersion tests using the testDispersion() function from {DHARMa}… first need to simulateResiduals
par(mfrow = c(1, 2))
P1M1_sim <- simulateResiduals(fittedModel = P1M1, plot = TRUE) #default is 250, authors may have increased because #s slightly dif
testDispersion(P1M1_sim)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.87988, p-value = 0.96
## alternative hypothesis: two.sided
NEED TO PULL OUT PVAL AND DEVIANCE RATIO FOR INTERACTION/INDEPENDENT DISPERSION TESTS Not getting exactly same value for deviance ratio… could have to do with how data was simulated?
check out this link
par(mfrow = c(1, 2))
P1M2_sim <- simulateResiduals(fittedModel = P1M2, plot = TRUE) #default is 250, authors may have increased because #s slightly dif
testDispersion(P1M2_sim)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.96419, p-value = 0.928
## alternative hypothesis: two.sided
Also got the same P val for this dispersion test but different deviance ratio (again prob because of simulation differences)
Prediction 2: Feeding & Travel
Lactating females don’t differ in feeding/travel time
‘We calculated our response variables by dividing the number of point samples that a given lactating female was engaged in feeding or travel, respectively, during each infant age class by the total number of good observations collected on that lactating female during that infant age class’
Feeding Model (P2M1)
first changed character values to factor since this is a diff file
#convert column 'id' from character to factor
feed_travel_social$ID <- as.factor(feed_travel_social$ID)
feed_travel_social$species <- as.factor(feed_travel_social$species)
feed_travel_social$age <- as.factor(feed_travel_social$age)
str(feed_travel_social)## 'data.frame': 34 obs. of 8 variables:
## $ species : Factor w/ 2 levels "bonobo","chimpanzee": 1 1 1 1 1 1 1 1 1 1 ...
## $ ID : Factor w/ 26 levels "BAH","Dju","DL",..: 2 11 13 15 17 24 25 26 19 22 ...
## $ age : Factor w/ 3 levels "0--18","18--36",..: 1 1 1 1 1 1 1 1 2 2 ...
## $ total : int 1618 3493 2794 2621 1425 2122 2717 1314 2196 1969 ...
## $ feed : int 662 1415 993 1089 383 707 1198 604 969 834 ...
## $ travel : int 316 386 471 370 202 400 586 225 490 478 ...
## $ social : int 272 720 579 367 95 367 366 140 279 223 ...
## $ social_adj: int 137 334 307 52 5 155 122 57 93 55 ...P2M1 <- glmmTMB(feed/total ~ species * age + (1 | ID), data = feed_travel_social, family = betabinomial(link = "logit"), weight = total)
P2M1_anova <- Anova(P2M1, type = c("III"), test.statistic = "Chisq")look at model/anova results and compare with paper
summary(P2M1)## Family: betabinomial ( logit )
## Formula: feed/total ~ species * age + (1 | ID)
## Data: feed_travel_social
## Weights: total
##
## AIC BIC logLik deviance df.resid
## 455.3 467.5 -219.7 439.3 26
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.03392 0.1842
## Number of obs: 34, groups: ID, 26
##
## Dispersion parameter for betabinomial family (): 56.6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.46398 0.11563 -4.013 6e-05 ***
## specieschimpanzee -0.11293 0.16074 -0.703 0.4823
## age18--36 0.11751 0.18741 0.627 0.5306
## age36--54 0.36437 0.17785 2.049 0.0405 *
## specieschimpanzee:age18--36 0.53117 0.28383 1.871 0.0613 .
## specieschimpanzee:age36--54 -0.05361 0.24885 -0.215 0.8294
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print(P2M1_anova)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: feed/total
## Chisq Df Pr(>Chisq)
## (Intercept) 16.1019 1 6.003e-05 ***
## species 0.4936 1 0.4823
## age 4.1999 2 0.1225
## species:age 4.3585 2 0.1131
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1pull out vals
P2M1_results <- tidy(P2M1_anova)
P2M1_x2 <- P2M1_results$statistic[4]
P2M1_df <- P2M1_results$df[4]
P2M1_p <- P2M1_results$p.value[4]Feeding Model (P2M2)
now looking at species and age as independent fixed effects
P2M2 <- glmmTMB(feed/total ~ species + age + (1 | ID), data = feed_travel_social, family = betabinomial(link = "logit"), weight = total)
P2M2_anova <- Anova(P2M2, type = c("II"), test.statistic = "Chisq")look at model/anova results and compare with paper
summary(P2M2)## Family: betabinomial ( logit )
## Formula: feed/total ~ species + age + (1 | ID)
## Data: feed_travel_social
## Weights: total
##
## AIC BIC logLik deviance df.resid
## 455.4 464.5 -221.7 443.4 28
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.028 0.1673
## Number of obs: 34, groups: ID, 26
##
## Dispersion parameter for betabinomial family (): 45.1
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.5071 0.1088 -4.662 3.12e-06 ***
## specieschimpanzee -0.0226 0.1256 -0.180 0.8572
## age18--36 0.3544 0.1512 2.344 0.0191 *
## age36--54 0.3189 0.1323 2.411 0.0159 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print(P2M2_anova)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: feed/total
## Chisq Df Pr(>Chisq)
## species 0.0324 1 0.85721
## age 8.3784 2 0.01516 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1age has significant effect on feeding model (supported by article)
pull out vals
P2M2_results <- tidy(P2M2_anova)
#pulling out values for species
P2M2_SPECIES_x2 <- P2M2_results$statistic[1]
P2M2_SPECIES_df <- P2M2_results$df[1]
P2M2_SPECIES_p <- P2M2_results$p.value[1]
#pulling out values for age
P2M2_AGE_x2 <- P2M2_results$statistic[2]
P2M2_AGE_df <- P2M2_results$df[2]
P2M2_AGE_p <- P2M2_results$p.value[2]###Traveling Model (P2M3)
Interact model first
P2M3 <- glmmTMB(travel/total ~ species * age + (1 | ID), data = feed_travel_social, family = betabinomial(link = "logit"), weight = total)
P2M3_anova <- Anova(P2M3, type = c("III"), test.statistic = "Chisq")look at model/anova results and compare with paper
summary(P2M3)## Family: betabinomial ( logit )
## Formula: travel/total ~ species * age + (1 | ID)
## Data: feed_travel_social
## Weights: total
##
## AIC BIC logLik deviance df.resid
## 391.5 403.7 -187.8 375.5 26
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.02284 0.1511
## Number of obs: 34, groups: ID, 26
##
## Dispersion parameter for betabinomial family (): 303
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.59316 0.07804 -20.415 <2e-16 ***
## specieschimpanzee -0.12538 0.11473 -1.093 0.274
## age18--36 0.14736 0.16062 0.917 0.359
## age36--54 0.15403 0.12064 1.277 0.202
## specieschimpanzee:age18--36 0.20379 0.25869 0.788 0.431
## specieschimpanzee:age36--54 -0.06929 0.16406 -0.422 0.673
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print(P2M3_anova)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: travel/total
## Chisq Df Pr(>Chisq)
## (Intercept) 416.7576 1 <2e-16 ***
## species 1.1943 1 0.2745
## age 2.6811 2 0.2617
## species:age 0.8496 2 0.6539
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1pull out vals
P2M3_results <- tidy(P2M3_anova)
P2M3_x2 <- P2M3_results$statistic[4]
P2M3_df <- P2M3_results$df[4]
P2M3_p <- P2M3_results$p.value[4]###Traveling Model (P2M4)
now looking at species and age as independent fixed effects
P2M4 <- glmmTMB(travel/total ~ species + age + (1 | ID), data = feed_travel_social, family = betabinomial(link = "logit"), weight = total)
P2M4_anova <- Anova(P2M4, type = c("II"), test.statistic = "Chisq")look at model/anova results and compare with paper
summary(P2M4)## Family: betabinomial ( logit )
## Formula: travel/total ~ species + age + (1 | ID)
## Data: feed_travel_social
## Weights: total
##
## AIC BIC logLik deviance df.resid
## 388.5 397.6 -188.2 376.5 28
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.01065 0.1032
## Number of obs: 34, groups: ID, 26
##
## Dispersion parameter for betabinomial family (): 203
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.60571 0.06880 -23.338 < 2e-16 ***
## specieschimpanzee -0.09255 0.08014 -1.155 0.24818
## age18--36 0.25062 0.09433 2.657 0.00789 **
## age36--54 0.10520 0.08471 1.242 0.21428
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print(P2M4_anova)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: travel/total
## Chisq Df Pr(>Chisq)
## species 1.3335 1 0.24818
## age 7.1529 2 0.02798 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1age has significant effect on feeding model (supported by article)
pull out vals
P2M4_results <- tidy(P2M4_anova)
#pulling out values for species
P2M4_SPECIES_x2 <- P2M4_results$statistic[1]
P2M4_SPECIES_df <- P2M4_results$df[1]
P2M4_SPECIES_p <- P2M4_results$p.value[1]
#pulling out values for age
P2M4_AGE_x2 <- P2M4_results$statistic[2]
P2M4_AGE_df <- P2M4_results$df[2]
P2M4_AGE_p <- P2M4_results$p.value[2]P2: Dispersion Tests
Next need to look at nonparametric dispersion tests using the testDispersion() function from {DHARMa}… first need to simulateResiduals
Feeding
par(mfrow = c(2, 2))
P2M1_sim <- simulateResiduals(fittedModel = P2M1, plot = TRUE) #feeding interaction
testDispersion(P2M1_sim)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.1414, p-value = 0.584
## alternative hypothesis: two.sidedP2M2_sim <- simulateResiduals(fittedModel = P2M2, plot = TRUE) #feeding independent
testDispersion(P2M2_sim)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.2964, p-value = 0.176
## alternative hypothesis: two.sided
NEED TO PULL OUT PVAL AND DEVIANCE RATIO FOR INTERACTION/INDEPENDENT DISPERSION TESTS Not getting exactly same value for deviance ratio… could have to do with how data was simulated?
Traveling
par(mfrow = c(1, 2))
P2M3_sim <- simulateResiduals(fittedModel = P2M3, plot = TRUE) #travel interaction
testDispersion(P2M3_sim)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.73323, p-value = 0.2
## alternative hypothesis: two.sidedP2M4_sim <- simulateResiduals(fittedModel = P2M4, plot = TRUE) #travel independent
testDispersion(P2M4_sim)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.80546, p-value = 0.368
## alternative hypothesis: two.sided
Also got the same P val for this dispersion test but different deviance ratio (again prob because of simulation differences)…
Prediction 3: Social Behavior
Lactating bonobos spend more time engaged in social interactions compared to chimpanzees
‘To test our third prediction that lactating bonobos spend more time engaged in social interactions, we ran two sets of models called Social Interactions and Adjusted Social Interactions. We calculated our response variable for social interactions by dividing the number of point samples that a given lactating female was engaged in social interactions during each infant age class by the total number of good observations collected on that lactating female during that infant age class. We calculated our response variable for adjusted social interactions by dividing the number of point samples that a given lactating female was engaged in social interactions with individuals other than her immature offspring during each infant age class by the total number of social interaction point samples collected on that lactating female during that infant age class.’
P3: Dispersion Tests
Next need to look at nonparametric dispersion tests using the testDispersion() function from {DHARMa}… first need to simulateResiduals
Social (P3M1 and P3M2)
par(mfrow = c(2, 2))
P3M1_sim <- simulateResiduals(fittedModel = P3M1, plot = TRUE) #social interaction
testDispersion(P3M1_sim)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.0657, p-value = 0.704
## alternative hypothesis: two.sidedP3M2_sim <- simulateResiduals(fittedModel = P3M2, plot = TRUE) #social independent
testDispersion(P3M2_sim)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.0331, p-value = 0.752
## alternative hypothesis: two.sided
NEED TO PULL OUT PVAL AND DEVIANCE RATIO FOR INTERACTION/INDEPENDENT DISPERSION TESTS Not getting exactly same value for deviance ratio… could have to do with how data was simulated?
Social Adjusted (P3M3 and P3M4)
par(mfrow = c(1, 2))
P3M3_sim <- simulateResiduals(fittedModel = P3M3, plot = TRUE) #travel interaction
testDispersion(P3M3_sim)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.90757, p-value = 0.968
## alternative hypothesis: two.sidedP3M4_sim <- simulateResiduals(fittedModel = P3M4, plot = TRUE) #travel independent
testDispersion(P3M4_sim)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.93123, p-value = 0.928
## alternative hypothesis: two.sided
Also got the same P val for this dispersion test but different deviance ratio (again prob because of simulation differences)…
Tables and Figures?
check out huxtable… glmmTMB not supported by stargazer or kable..
Social Interaction (P3M1)
first test interaction of species and age when looking at social behav
look at model/anova results and compare with paper
pull out vals